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We present results for the isovector and flavor diagonal tensor charges glf~'^, gif, 5ti s^nd 
needed to probe novel tensor interactions at the TeV scale in neutron and nuclear /3-decays and 
the contribution of the quark electric dipole moment (EDM) to the neutron EDM. The lattice 
QCD calculations were done using nine ensembles of gauge configurations generated by the MILC 
collaboration using the HISQ action with 2-I-1-I-1 dynamical flavors. These ensembles span three 
lattice spacings a ~ 0.06, 0.09 and 0.12 fm and three quark masses corresponding to the pion masses 
Mt, « 130, 220 and 310 MeV. Using estimates from these ensembles, we quantify all systematic 
uncertainties and perform a simultaneous extrapolation in the lattice spacing, volume and light quark 
masses for the connected contributions. The final estimates of the connected nucleon (proton) tensor 
charge for the isovector combination is = 1.020(76) in the MS scheme at 2 GeV. The additional 
disconnected quark loop contributions needed for the flavor-diagonal matrix elements are calculated 
using a stochastic estimator employing the truncated solver method with the all-mode-averaging 
technique. We find that the size of the disconnected contribution is smaller than the statistical error 
in the connected contribution. This allows us to bound the disconnected contribution and include 
it as an additional uncertainty in the flavor-diagonal charges. After a continuum extrapolation, 
we find gf — 0.774(66), gtf = —0.233(28) and = 0.541(67). The strangeness tensor charge, 
that can make a significant contribution to the neutron EDM due to the large ratio mslmu,d, is 
gx = 0.008(9) in the continuum limit. 
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I. INTRODUCTION 

Precise estimates of the matrix elements of the 
isoscalar and isovector tensor bilinear quark operators 
are needed to obtain bounds on new physics from preci¬ 
sion measurements of /3-decays and limits on the neutron 
electric dipole moment (nEDM). The isovector charge, 
is needed to probe novel tensor interactions in the 
helicity-flip part of the neutron decay distribution [T] 
while the isoscalar charges are needed to quantify the 
contribution of the quark EDM to the nEDM and set 
bounds on new sources of CP violation. In this paper, 
we give details of the simulations of lattice QCD using 
the clover-on-HISQ approach to provide first principle 
estimates of these matrix elements with control over all 
sources of systematic errors. 

Lattice QCD analysis of isovector charges of nucleons 
is well-developed (See the recent reviews [2H1])- In this 
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work we present precise estimates of this dominant con¬ 
tribution, given by the connected diagrams, to the tensor 
charges, i.e., the insertion of the zero-momentum tensor 
operator in one of the three valence quarks in the nucleon. 
Calculation of the isoscalar charges is similar except that 
it gets additional contributions from contractions of the 
operator as a vacuum quark loop. This is called the dis¬ 
connected contribution as the quark loop and nucleon 
propagator interact only through the exchange of gluons. 
The statistical signal in the disconnected term is weak, 
so it is computationally much more expensive. We find, 
on the four ensembles analyzed, that the disconnected 
contributions of light quarks is small and in most cases 
are consistent with zero within errors. We, therefore, use 
the largest of these estimates to bound the disconnected 
contribution and include it as a systematic uncertainty in 
the presentation of the final results. Similarly, using five 
ensembles we show that the disconnected contribution of 
the strange quark, also needed for the nEDM analysis, 
is even smaller but we are able to extract a continuum 
limit estimate [S]. 

Throughout the paper, we present results for the ten¬ 
sor charges of the proton, which by convention are called 
nucleon tensor charges in literature. Results for the neu¬ 
tron are obtained by the u ^ d interchange. This paper 
is organized as follows. In Section |llj we describe the 
parameters of the gauge ensembles analyzed, the lattice 
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methodology, fits used to extract matrix elements within 
the ground state of the nucleon and the renormalization 
of the operators. We discuss the calculation of the con¬ 
nected diagram in section [ml and of the disconnected 
contribution in Section IV Our final results are presented 
in Section and we end with conclusions in Section IVII 
In the Appendix we present a summary of the control 
over systematics of existing lattice calculations using the 
FLAG quality criteria [H]. 


II. LATTICE PARAMETERS AND 
METHODOLOGY 

In this section we provide an overview of the calcula- 
tional details. These include a description of the gauge 
ensembles analyzed, a short review of the operators used 
to calculate the two-point and three-point correlation 
functions using clover fermions, the fit ansatz used to 
extract the desired matrix elements from the correlation 
functions and estimates of renormalization constants us¬ 
ing the RI-sMOM scheme. 


A. Lattice Parameters 

In order to obtain estimates with a desired precision, 
it is important to quantify all sources of systematic er¬ 
rors. For matrix elements between nucleon ground states, 
these include excited state contamination, finite lattice 
volume, operator renormalization, discretization effects 
at finite lattice spacing and extrapolations from heavier 
u and d quarks. Since lattice generation is very expen¬ 
sive, it was, therefore, expedient to use a set of existing 
gauge ensembles that cover a sufficiently large range in 
lattice spacing and light quark mass to study the con¬ 
tinuum and chiral behavior. The only set available to us 
that meets our requirements are the ensembles generated 
using Nf = 2 -I-1 -I-1 flavors of highly improved staggered 
quarks (HISQ) [7] by the MILC collaboration [5]. The 
parameters of the nine ensembles used in this study are 
given in Table |l] In this paper we show that these ensem¬ 
bles allow us to address issues of statistics, excited state 
contamination, lattice volume, lattice spacing and the 
chiral behavior in the calculation of the tensor charges. 

Staggered fermions have the advantage of being com¬ 
putationally cheaper and preserve an important remnant 
of the continuum chiral symmetry. Their disadvantage 
is that the spectrum has a four-fold doubling in the con¬ 
tinuum limit. This doubling symmetry (called the taste 
symmetry) is broken at finite lattice spacing and this 
breaking introduces additional lattice artifacts. Due to 
taste mixing, staggered baryon interpolating operators 
couple, in general, to a combination of octet (the nu¬ 
cleon) and the decuplet (Delta) states. Furthermore, 
these baryon operators couple to both parity states in 
addition to all radial excitations of these. Thus, baryon 
correlation functions are more complicated to analyze 


compared to Wilson-type fermions, as they have a weaker 
statistical signal, the consequences of taste mixing has to 
be resolved and one has to take into account the oscillat¬ 
ing signal due to contributions from both parity states. 
Since having a good statistical signal is a prerequisite 
to quantifying the various sources of systematic errors, 
we have chosen to construct correlation functions using 
Wilson-clover fermions, as these preserve the continuum 
spin structure. This mixed-action, clover-on-HISQ, ap¬ 
proach, however, leads to a nonunitary lattice formula¬ 
tion and at small, but a priori unknown, quark masses 
suffers from the problem of exceptional configurations 
discussed next. 

Exceptional configurations are ones in which the clover 
Dirac operator evaluated on HISQ configurations has 
near zero modes. As a result, on such configurations 
the inversion of the clover Dirac operator, which gives 
the Feynman propagator, can fail to converge and/or the 
corresponding correlation functions have an exception¬ 
ally large amplitude depending on the proximity to the 
zero mode. The presence of exceptional configurations 
biases the results or gives rise to unphysically large fluctu¬ 
ations and invalidates the results. In any lattice analysis 
based on a unitary formulation, such as HISQ-on-HISQ 
or clover-on-clover, such configurations are suppressed in 
the lattice generation. Given an appropriately generated 
ensemble of HISQ configurations, there is no basis for ex¬ 
cluding any configuration from the ensemble average in 
a clover-on-HISQ calculation. Thus, these calculations 
should be done only on ensembles without any excep¬ 
tional configurations. 

The presence of such exceptional configurations in a 
clover-on-HISQ analysis is expected to increase on de¬ 
creasing the quark mass at fixed lattice spacing and in¬ 
crease with the lattice spacing at fixed quark mass, i.e., 
the coarser the configurations, the more likely they are. 
Consequently, smearing techniques used to reduce short 
distance lattice artifacts also reduce the probability of en¬ 
countering exceptional configurations. To reduce lattice 
artifacts, we applied HYP smearing [9] to all HISQ con¬ 
figurations. To check for exceptional configurations on 
these HYP smeared lattices, we monitor the convergence 
of the quark propagator and the size of fluctuations in 
correlations functions on each configuration. These tests 
provided evidence for exceptional configurations on the 
a = 0.15 fm ensembles and on the a = 0.12 fm ensem¬ 
ble with Mtt « 130 MeV. Consequently, these ensembles, 
also generated by the MILC collaboration, were excluded 
from our analysis. An earlier discussion regarding ex¬ 
ceptional configurations on these ensembles is given in 
Ref. [TU]. To reiterate, the nine ensembles used in this 
study and described in Table |l] did not present evidence 
of an exceptional configuration. 

The parameters used in the analysis with clover 
fermions are given in Table [H] The Sheikholeslami- 
Wohlert coefficient used in the clover action is fixed to 
its tree-level value with tadpole improvement, = 1 /wq 
where uq is the tadpole factor of the HYP smeared HISQ 
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Ensemble ID 

a (fm) 

(MeV) 

xT 


al2m310 

A 

0.1207(11) 

305.3(4) 

24® X 64 

4.55 

al2m220S 

< 

0.1202(12) 

218.1(4) 

24® X 64 

3.29 

al2m220 

0 

0.1184(10) 

216.9(2) 

32® X 64 

4.38 

al2m220L 

V 

0.1189(09) 

217.0(2) 

40® X 64 

5.49 

a09m310 

<0> 

0.0888(08) 

312.7(6) 

32® X 96 

4.51 

a09m220 

☆ 

0.0872(07) 

220.3(2) 

48® X 96 

4.79 

a09ml30 

0 

0.0871(06) 

128.2(1) 

64® X 96 

3.90 

a06m310 

□ 

0.0582(04) 

319.3(5) 

48® X 144 

4.52 

a06m220 

0 

0.0578(04) 

229.2(4) 

64® X 144 

4.41 


TABLE I. Parameters of the (2+1+1) flavor HISQ lattices 
generated by the MILC collaboration and analyzed in this 
study are quoted from Ref. [8]. Symbols used in plots are 
defined along with the ensemble ID. Finite size effects are an¬ 
alyzed in terms of M^L with the clover-on-HISQ defined 
in Table HI 


lattices. 

The masses of light clover quarks were tuned so that 
the clover-on-HISQ pion masses (see Table [ll|) match 
the HISQ-on-HISQ Goldstone ones, given in Ta¬ 

ble The strange quark mass nis is tuned so that the 

resulting clover-on-HISQ Mgg = , where 

w!l/m\ is the ratio of bare strange and light quark masses 
used in the HISQ generation, and is 5 for m3I0 lattices, 
10 for m220 lattices and 27 for ml30 lattices. The result¬ 
ing estimates for mg are consistent with those obtained 
by matching to the kaon mass m- 

All fits in to study the chiral behavior are made us¬ 
ing the clover-on-HISQ as correlation functions and 
thus the observables have a greater sensitivity to it. Per¬ 
forming hts using the HISQ-on-HISQ values of did 

not change the estimates significantly. 

Estimates of nucleon charges and form-factors based on 
lower statistics subset of data presented here for the two 
ensembles al2m310 and al2m220 have been published 
in [To]. Results for the tensor charges presented in this 
paper supersede those earlier estimates. 


B. Lattice Methodology 

The two-point and three-point nucleon correlation 
functions at zero momentum are defined as 

X 

x,x' 

where a and 13 are the spinor indices. The source time 
slice is translated to to = 0, t is the sink time slice, and r 
is the time slice at which the bilinear operator Op(a;) = 
q{x)Tq(x) is inserted. The Dirac matrix T is 1, 74 , 7^75 


and 7 i 7 j for scalar (S), vector (V), axial (A) and tensor 
(T) operators, respectively. In this paper, subscripts i 
and j on gamma matrices run over {1,2,3}, with i < j. 
The interpolating operator used to create/annihilate the 
nucleon state, %, is 


X{x) = 


^abc 


(x)C'75-(l±74)<72(a;) 




( 3 ) 


with color indices {a, b, c}, charge conjugation matrix C, 
and qi and 52 denoting the two different flavors of light 
quarks. The nonrelativistic projection (1 ± 74)/2 is in¬ 
serted to improve the signal, with the plus and minus 
sign applied to the forward {t > 0) and backward {t < 0) 
propagation, respectively. 

The nucleon charges gf are dehned as 

{N{p,s)\0^\N{p,s)) = gl,Us{p)rus{p) (4) 
with spinors satisfying 

Us {p)us (p) = |< + rriN . (5) 

S 

These charges, g'f, can be extracted from the ratio of the 
projected three-point function to the two-point function 
for f > T > 0 


Rr{t,T) 


^ (Tr[PrCgPyf,T)]) 

(Tr[P2ptC2pHt)]) 

—^ Tr [Pr(l + 74)r(l + 74 )] gf. 


( 6 ) 


Here V 2 pt = (l + 74)/2 is used to project out the positive 
parity contribution and T^r is dehned below. Note that 
the ratio in Eq. ®> becomes zero if T anti-commutes with 
74 , so only r = 1, 74 , 7 i 75 and 7 ^ 7 ^ can survive. In this 
paper, we present results for only the tensor channel, 
for which we can demonstrate control over all systematic 
errors. 

On inserting a bilinear quark operator between the nu¬ 
cleon states to construct the three-point function, one 
gets two classes of diagrams: (i) the bilinear operator is 
contracted with one of the three valence quarks in the 
nucleon, as shown in the left diagram of Fig. [T and (ii) 
the bilinear operator is contracted into a quark loop that 
is correlated with the nucleon two-point function through 
the exchange of gluons, as shown in the right diagram of 
Fig.[T| These are called the connected and disconnected 
diagrams, respectively. 

The disconnected part of Eq. can be written as 
= (y^Tr[M"i(r,x;T,x)r]) 

X 

{Tr[VrC^^\t)] Ex Tr[M-i(r, x; r, x)r]) 

(Tr[7^2ptC2pt(t)]) ’ ^ > 

where M is the Dirac operator. Note that the first term 
of the right-hand side is zero when E 7 ^ 1 , so it does not 
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ID 

mi 

ms 

csw 

Mr' (MeV) 

Smearing 

al2m310 

-0.0695 

-0.018718 

1.05094 

310.2(2.8) 

{5.5, 70} [{5.5, 70}] 

al2m220S 

-0.075 

— 

1.05091 

225.0(2.3) 

{5.5, 70} 

al2m220 

-0.075 

-0.02118 

1.05091 

227.9(1.9) 

{5.5, 70} [{5.5, 70}] 

al2m220L 

-0.075 

— 

1.05091 

227.6(1.7) 

{5.5, 70} 

a09m310 

-0.05138 

-0.016075 

1.04243 

313.0(2.8) 

{5.5, 70} [{6.0, 80}] 

a09m220 

-0.0554 

-0.01761 

1.04239 

225.9(1.8) 

{5.5, 70} [{6.0, 80}] 

a09ml30 

-0.058 

— 

1.04239 

138.1(1.0) 

{5.5, 70} 

a06m310 

-0.0398 

-0.01841 

1.03493 

319.6(2.2) 

{6.5, 70} [{6.5, 80}] 

a06m220 

-0.04222 

— 

1.03493 

235.2(1.7) 

{5.5, 70} 


TABLE II. The parameters used in the calculation of clover propagators. The hopping parameter k in the clover action is 
given by + 4). rris is needed for the calculation of the strange quark disconnected loop diagram. The Gaussian 

smearing parameters are defined by {cr, A^kg} where Nkg is the number of applications of the Klein-Gordon operator and 
the width is controlled by the coefficient cr, in Ghroma convention. Smearing parameters used in the study of disconnected 
diagrams are given within square-brackets, mi is tuned to achieve and mg is tuned so that Mgg = y/ 

The error in the pion mass M,r is governed mainly by the uncertainty in the lattice scale given in Table [I] 


contribute to the tensor charges. High precision mea¬ 
surements of Eq. Q requires improving the signal in the 
second term, i.e., the correlation between the nucleon 
two-point functions and the quark loop Tr[M“^r] 
as discussed in Sections llV Al and IIV Bl 

The charges gf are extracted from the ratio Rr by 
appropriate choice of the projection operator Vr- For 
the calculation of connected contribution we use a single 
projection operator Vr — ’P 2 pt(l + * 7573 ) to extract all 
four tensor structures at the same time as the projection 
is done at the time of the calculation of the sequential 
propagator. For the disconnected diagram, the projec¬ 
tion operator is part of the final trace with the two-point 
function, so there is no additional cost to using tensor 
specific projectors. We, therefore, use 

Vi = 7^2pt, 

^7i75 = ^2pt757i (* = 1, 2, 3), 

"^7172 = ^2pt7573, (8) 

^7273 '^2pt757l) 

'^7173 “ '^2pt7275j 

which make Rr{t, t) — >• gr for f 3> t 0 as they satisfy 
iTr[Pr(l+74)r(l+74)] = 1. (9) 

For the disconnected diagrams, the statistics in the 
two-point function are improved by including the back¬ 
ward propagating baryons [t < 0) from each source point. 
In this case V^^*" = (1 — 74)/2 is used to project out the 
negative parity state and we multiply the tensor projec¬ 
tion operators V-y-^- in Eq. ([^ by (—1) to match the 
convention used for forward propagation. 

The calculation of the two-point and the connected 
three-point functions is carried out using the method de¬ 
scribed in our previous study. Ref. [10]. These calcu¬ 
lations use the Chroma software package [l2|. Part of 
the calculations are done on clusters with graphic pro¬ 
cessing units (GPUs) using QUDA library j^. The 


source and sink baryon operators are constructed using 
smeared quark propagators to reduce the contamination 
from the excited states. We use gauge-invariant Gaus¬ 
sian smeared sources to improve the overlap with the 
ground state. Smearing is done by applying the three- 
dimensional Klein-Gordon operator a fixed number of 
times (1 —(t^V^/(4A^kg))^’^‘^- The smearing parameters 
{cr, IVkg} for each ensemble are given in Table [Hj 

To calculate the connected three-point function, we an¬ 
alyze configurations in sets of four measurements, i.e. we 
generate four independent propagators on each con¬ 
figuration using smeared sources on four maximally sep¬ 
arated time slices For each S°, the same smear¬ 

ing operation is applied at all time slices to create the 
smeared sink, and the two-point correlation function is 
calculated using these smeared-smeared propagators Sss- 
Each of these four smeared propagators are used to con¬ 
struct sources for u and d quark propagators with the 
insertion of zero-momentum nucleon state at time slices 
displaced by a fixed Gep from the four source time slices 
These u and d sequential sources (generated sep¬ 
arately) at + tsep are smeared again. The final 

coherent sequential propagator is then calculated 
using the sum of these four smeared sources, be., the se¬ 
quential propagator from the four time slices is calculated 
at one go. The connected three-point functions, over the 
four regions to + tsep, are then constructed 

by inserting the bilinear operator between each of the 
original individual propagator S° from and the co¬ 

herent sequential propagator from '’"+tsep M- 
The assumption one makes by adding the four sources to 
produce a single coherent sequential propagator is that 
the entire contribution to the three-point function in any 
one of the four intervals between and +tsep is 

from baryon insertion at -I- tgep and the contribution 
of baryon sources at the other three time slices + tgep 
in S'®®'* goes to zero on averaging over the gauge configu¬ 
rations. The coherent sequential source method has the 
advantage that the insertion of operators with different 
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tensor structures and various momenta and for all four 
source positions can be done at the same time with tiny 
computational overhead. 

To study and quantify the excited state contamination, 
we repeat the calculation for multiple source-sink separa¬ 
tions, tsep- Separate sequential u and d propagators are 
calculated for each 4ep analyzed. Thus, the total number 
of inversions of the Dirac operator are 4 -|- 2 x for 
each set of four measurements on each configuration. Our 
choices of tsep and the number of measurements made 
on each ensemble (number of configurations times the 
number of sources on each configuration) are given in 

Table mu 

The calculation of disconnected quark loop diagrams 
using stochastic methods have a poor signal and requires 
very high statistics. Because of the computational cost, 
the calculations with light quarks have been done on the 
three heaviest, « 310 MeV, and the al2m220 en¬ 
sembles; and on five ensembles for the strange quark as 
listed in Table IVIIII For the evaluation of the discon¬ 
nected diagrams, we obtain a stochastic estimate using 
the truncated solver method (TSM) [lU [16] with the 
all-mode-averaging (AMA) technique [T7] as described 
in Section HVl 


C. Fits to Correlation Functions 


To extract the desired nucleon charges, the matrix el¬ 
ements of the bilinear quark operators need to be cal¬ 
culated between ground state nucleons. On the lattice, 
however, any zero-momentum correlation function de¬ 
fined in Eq. ([^ using the nucleon interpolation operator 
defined in ([^, has a coupling to the ground state nucleon, 
all radially excited states, and multiparticle states with 
the same quantum numbers. Operators constructed us¬ 
ing appropriately tuned smeared sources reduce the cou¬ 
pling to excited states but do not eliminate it. We discuss 
two synergistic strategies for removing the remaining ex¬ 
cited state contamination based on the fact that in Eu¬ 
clidean time, the contributions from the excited states 
are exponentially suppressed as (i) the distance between 
the source/sink and the inserted operator increases and 
(ii) as the mass gap, Mexcited — Afg, increases. First, 
one can increase tsep- Unfortunately, the signal also de¬ 
creases exponentially as tsep is increased so one is forced 
to compromise. In fact, the values of tsep we have used re¬ 
flect this compromise based on the anticipated statistics 
on each ensemble. Additionally, one can include excited 
states in the analysis of Eq. © as discussed below. 

We include one excited state in the analysis of the two- 
and three-point functions. For operator insertion at zero 


momentum, the data are fit using the ansatze 
C^p\tf,U) = 

^ (iq) 

|AoP(0|Or|0)e-^A(t/-q) + 


where the source positions are shifted to = 0 and 
tf = tsep. The states |0) and |1) represent the ground 
and “first” excited nucleon states, respectively. The four 
parameters, Mq, Mi, Ap and Ai are estimated first from 
the two-point function data. We find that the extraction 
of Mg and Ao is stable under change of the fit range, 
while that of Mi and Ai is not. We, therefore, choose 
the largest range, requiring that the values of, and the 
errors in, all four parameters do not jump by a large 
amount on changing the fit range. In all these fits, we 
find Ml « 2Mq, so it should be considered an effective 
excited state mass as it is much larger than the 7V(1440) 
excitation. The results of these best fits are given in Ta¬ 
ble HVl 

We performed two independent measurements of the 
two-point functions, and the corresponding Mg, Mi, Aq 
and Ai are given in Table |TVj The second set of measure¬ 
ments were obtained during the calculation of the discon¬ 
nected diagrams using the AMA error reduction method 
discussed in Section lTV AI We find that the two estimates 
are consistent within errors indicating no remaining bias 
with our choice of parameters for the AMA. 

Fits using the ansatz for the three-point function given 
in Eq. (II) are used to isolate the two unwanted matrix 
elements(0|Or|l) and (l|Or|l)- We find that the mag¬ 
nitude of (0|(!lr|l) is about 16% of (0|C>r|0) and is deter¬ 
mined with about 20% uncertainty on all the ensembles, 
whereas |(l|Or|l)| (0|Or|0), but has 0(100%) errors. 

Ideally, equally precise data should be generated at each 
value of tsep- In our analysis, however, the same num¬ 
ber of measurements have been made for all tsep on each 
ensemble, so errors increase with tgep ns shown in Fig. [21 

We reduce the contamination from higher excited 
states in these fits to the two- and three-point functions 
by excluding data points overlapping with, and adjacent 
to, the source and sink time slices at which the excited 
state contamination is the largest. For uniformity, we 
exclude 2, 3, 4 time slices on either end of the inter¬ 
val tsrc to tsrc + tsep In the three-point function for the 
a = 0.12, a = 0.09 and a = 0.06 fm ensembles, respec¬ 
tively. In physical units, these excluded regions corre¬ 
spond to roughly the same distance. This range of ex¬ 
cluded points is consistent with the starting time slice of 
the fits to the two-point correlators given in Table |IV[ 
be., the time beyond which a two-state fit captures the 
two-point function data. Changing the exclusion time 
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FIG. 1. The connected (left) and disconnected (right) three-point diagrams needed to calculate the matrix elements of bilinear 
quark operators in the nucleon state. 


ID 

^sep / Ct 

A^conf 

^meas 

al2m310 

{8,9,10,11,12} 

1013 

8104 

al2m220S 

{8,10,12} 

1000 

24000 

al2m220 

{8,10,12} 

958 

7664 

al2m220L 

10 

1010 

8080 

a09m310 

{10,12,14} 

881 

7048 

a09m220 

{10,12,14} 

890 

7120 

a09ml30 

{10,12,14} 

883 

7064 

a06m310 

{16,20,22,24} 

1000 

8000 

a06m220 

{16,20,22,24} 

650 

2600 


TABLE III. The values of source-sink time separations 
(taep/a) used, the total number of configurations analyzed 
(A^conf) and measurements made (Ameas) for the two-point 
and connected three-point function calculations. 


slice values to 3, 4, and 6 in both the fits changed the 
final estimates of the charges by less than Itr. 

Including a second exited state in the analysis would 
increase the number of matrix elements to be estimated 
from the three-point function by three. Given that their 
contribution would be smaller still, much higher statistics 
than generated for this study would be needed. This is 
confirmed in practice; the data are well fit with just the 
one excited state ansatz and there is no sensitivity left to 
resolve three additional small parameters. This analysis, 
using the ansatze given in Eq. 0 and fitting the data at 
all tsep simultaneously, was called the two-simRR method 
in Ref. [TO] . 

Our overall conclusion is that, using the values of tgep 
and the statistics for each ensemble given in Table m 
and assuming that only one excited state gives a signifi¬ 
cant contribution, we are able to isolate and remove this 
contamination as illustrated in Fig. It turns out that 
on all nine ensembles, the excited state contamination 
for gx is small. It is worth remarking that this is not the 
case for gA as discussed in m- 


D. Renormalization of Operators 

The calculation of the renormalization constants Zp of 
the quark bilinear operators in the RI-sMOM scheme [El 
[E] has been done on five ensembles: al2m310, al3m220, 
a09m310, a09m220 and a06m310. In order to translate 
the lattice results to the continuum MS scheme at a fixed 
scale, say g, = 2 GeV, used by phenomenologists we fol¬ 
low the procedure described in Ref. m- To summarize, 
the RI-sMOM estimate obtained at a given lattice four- 
momentum is first converted to the MS scheme at 
the same scale (horizontal matching) using the one-loop 
perturbative matching. This value is then run in the con¬ 
tinuum in the MS scheme to the fixed scale, 2 GeV, using 
the two-loop anomalous dimension. 

Ideally, one would like to establish a window A ^ g <C 
c/a in the RI-sMOM scheme in which the Zp scales ac¬ 
cording to perturbation theory. Here A is an infrared 
scale below which nonperturbative effects are large and 
c/a represents the cutoff scale beyond which lattice dis¬ 
cretization effects are large. The value of c is o priori 
unknown and the expectation is that it is 0(1). Within 
this window, the scaling of Zt with gets contributions 
from both the anomalous dimension of the operator and 
the running of the strong coupling constant Og. If this 
scaling is consistent with that predicted by perturbation 
theory, then estimates within this window would con¬ 
verge to a constant value independent of after conver¬ 
sion to MS scheme and run to a fixed scale, 2 GeV. As 
discussed in [E] , smearing the lattice to reduce the ultra¬ 
violet noise in the measurements also reduces the upper 
cutoff c/a for the calculation of the renormalization con¬ 
stants, and a priori, we again do not know by how much 
smearing shrinks the desired window or whether it totally 
eliminates it on the various 0.06—0.12 fm lattices we have 
analyzed. Below we summarize the tests performed and 
state the results. 

• We first test the data for the Z's in the RI-sMOM 
scheme to see if they exhibit the desired perturba- 
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FIG. 2. The data for g^~ and the results of the simultaneous fit using multiple tsep using the ansatz given in Eq. (111 to isolate 
the excited state contribution. The seven figures are arranged as follows: the « 310 MeV ensembles (top), « 220 MeV 
ensembles (middle) and the ~ 130 MeV ensemble (bottom). The solid black line and the grey band are the ground state 
(taep —>■ oo) estimate and error. The fits evaluated for different taep are also shown by solid lines. 


live behavior for HYP smeared lattices by calculat¬ 
ing the logarithmic derivative of Z{q^) and compar¬ 
ing it to the anomalous dimension. The data show 
evidence of such a window in the calculation of the 
vector, axial and tensor renormalization constants, 
but not for the scalar. In this paper, we only need 
Zt and Zy, so we next describe how we obtained 
final estimates for these and assigned a conservative 
error that covers the various sources of systematic 
uncertainties. 

• We find that the ratios of renormalization con¬ 
stants, ZylZy^ have less fluctuations and are flat¬ 
ter in as illustrated in Fig. This improvement 
is presumably due to the cancellation of some of 
the systematic uncertainties in the ratio, includ¬ 
ing, for example, those due to the breaking of the 
continuum Lorentz symmetry to the hypercubic ro¬ 
tation group on the lattice that impacts the calcu¬ 


lation of the Z-factors. On each ensemble, the fi¬ 
nal renormalized charges can be constructed from 
these ratios as {ZyjZv) x (ST/Sy””^) using the iden¬ 
tity Zvgy~‘^ = 1. Because of the better signal and 
resulting fits, we use the estimates from the ratios 
method for our central values and include the dif¬ 
ference between these and estimates from the direct 
calculation, Z^gr, as an estimate of the systematic 
error. 

• To take into account the remaining dependence on 
q^ of the estimates in MS scheme at 2 GeV, we carry 
out the two analysis strategies proposed in Ref. [10] . 
In the first, we obtain the value and error from 
the fit to the data using the ansatz c/q"^ + Z + aq. 
We find that these fits capture the data and the 
extrapolations Z -\- aq are shown as dashed lines in 
Fig.i 

• A slightly modified version of the second method 





































































ID 

Fit Range 

aMo 

aMi 

Ao X 10“ 

Ai X 10“ 

Fit Range 

aMo 

aMi 

Ao X 10“ 

Ai X 10“ 

al2m310 

al2m220S 

2-15 

2-15 

0.6669(53) 

0.6233(55) 

1.36(11) 

1.42(13) 

6.57(27) 

6.58(26) 

6.28(61) 

6.94(93) 

2-15 

0.6701(16) 

1.471(45) 

6.845(82) 

6.88(35) 

al2m220 

al2m220L 

2-15 

2-15 

0.6232(49) 

0.6046(71) 

1.45(15) 

1.16(12) 

6.58(24) 

5.68(37) 

6.8(1.1) 
5.63(51) 

2-15 

0.6124(17) 

1.294(37) 

6.070(91) 

6.34(23) 

a09m310 

3-20 

0.4965(46) 

0.938(57) 

14.12(75) 

17.4(1.1) 

3-20 

0.4973(12) 

0.971(22) 

2.215(31) 

2.374(74) 

a09m220 

a09ml30 

3-20 

3-20 

0.4554(45) 

0.4186(76) 

0.925(53) 

0.834(61) 

12.13(61) 

9.74(89) 

18.5(1.3) 

17.2(1.0) 

3-20 

0.4524(24) 

0.877(34) 

1.812(56) 

2.29(10) 

a06m310 

a06m220 

4- 30 

5- 30 

0.3245(30) 

0.3166(66) 

0.617(18) 

0.644(54) 

0.566(30) 

13.0(1.5) 

1.439(42) 

38.5(5.4) 

4-30 

0.3283(15) 

0.630(10) 

0.609(15) 

1.513(29) 


TABLE IV. The nucleon ground and first excited state masses and the corresponding amplitudes obtained from a two-state 
fit to the nucleon two-point correlation function on each ensemble. The second set of estimates on the right are from an 
independent calculation performed to calculate the disconnected diagrams using the AMA with 64 LP measurements (96 LP 
for a06m310). All errors are estimated using the single-elimination Jackknife method using uncorrelated fits. 


is used: we now choose the in the RI-sMOM 
scheme by the condition qia — s\n{qia) = 0.05 based 
on bounding the discretization error and the error 
in Z is estimated from the spread in the data over 
a range in q^ about this point. This choice corre¬ 
sponds to = 5, 9 and 21 GeV^ for the a = 0.12, 
0.09 and 0.06 fm ensembles, respectively. The cor¬ 
responding ranges for determining the error were 
taken to be 4-6, 8-10 and 18-24 GeV^ over which 
the data show a reasonably flat behavior as shown 
in Fig. ^ 

• For the final estimates we take the average of the 
two methods. The error is taken to be half the 
difference, and rounded up to be conservative. 

• On the ensembles at the two lattice spacings a = 
0.12 and 0.09 fm, we found no significant difference 
in the estimates of the renormalization constants 
for the two different quark masses (Mtt = 310 and 
220 MeV ensembles). A common fit captured both 
data sets, as shown in Fig. and was used to ex¬ 
tract our “quark mass independent” estimates. 

• The entire calculation, matching to the MS scheme, 
running to 2 GeV and the final fits for the two 
strategies, was done using 200 bootstrap samples 
because the number of configurations analyzed in 
ensembles al2m310 and al2m220 {aOQmSlO and 
a09m220) are different. 

The final mass-independent renormalization constants 
at the three lattice spacings needed to construct the 
renormalized charges in the two ways: (i) .^r5r and (ii) 
from the product of the ratios {Zy-/Zv) x with 

the identity Zvgy~‘^ = 1 are given in Table |v| The errors 
in the renormalization factors, Zt and Zt/^, are added 
in quadratures to those in the extraction of the bare nu¬ 
cleon charges and respectively, to get 

the final estimates of the renormalized charges given in 
Table rvin 


E. Statistical analysis of two-point and three-point 
functions 

We carried out the following statistical analyses of the 
data on each ensemble to look for anomalies. We divided 
the data for a given ensemble into bins of about 1000 mea¬ 
surements (by source points and by configuration gen¬ 
eration order) to test whether the ensembles consist of 
enough independent configurations. For bin sizes > 5 
configurations, the errors in the mean decreased as y/N, 
i.e., consistent with our analysis of the autocorrelation 
coefficient of about 5 configurations (about 25 molecu¬ 
lar dynamics steps). Also, the error computed with data 
averaged over S source points on each configuration is 
smaller by y/S compared to the error in the data from 
any one of the source position. 

Estimates from bins of about 1000 measurements, how¬ 
ever, fluctuated by up to 3tT in some cases. This variation 
is much larger than expected based on the bin sizes. To 
determine whether the data in the various bins satisfy 
the condition of being drawn from the same distribu¬ 
tion, we performed the Kolmogorov-Smirnov (K-S) test 
on quantities that have reasonable estimates configura¬ 
tion by configuration, for example, the isovector vector 
charge gy~'^ and the value of the two-point function at a 
given time separation. The K-S test showed acceptable 
probability of the various bins being drawn from the same 
distribution. Histograms of the data showed no long tails 
in the distribution but exhibit variations in the sample 
distribution that becomes increasingly Gaussian as the 
bin size was increased to the full sample size. 

We find these 2 — 3(7 fluctuations both when the data 
are binned by the source position and when the con¬ 
figurations are divided into two halves according to the 
molecular dynamics generation order. Such fluctuations 
are apparent in the aOSmSl 0 and the a06m220 ensemble 
data. Comparing the data for different charges (axial, 
scalar and tensor), we found that the effect is least sig¬ 
nificant (less than Icr) for the tensor charge and worst 
for the vector charge gy] it is, presumably, most evident 
in gy because it has the smallest statistical errors. We 
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ID 

Zt 

Zv 

Zt jZy 

Zt 

Zv 

ZtIZv 

Zt 

Zv 

Zt jZv 

al2 

0.898(4) 

0.890(4) 

1.003(3) 

0.995(10) 

0.918(12) 

1.073(2) 

0.95(5) 

0.90(2) 

1.04(4) 

a09 

0.962(6) 

0.911(9) 

1.045(5) 

1.026(6) 

0.938(6) 

1.089(2) 

0.99(4) 

0.925(15) 

1.07(3) 

a06 

1.005(6) 

0.931(4) 

1.072(5) 

1.071(5) 

0.961(5) 

1.1134(6) 

1.04(4) 

0.945(15) 

1.09(3) 


TABLE V. The mass independent renormalization constants Zt, Zy and the ratio ZtIZv in the MS scheme at 2GeV at the 
three values of the lattice spacings used in onr calculations. These estimates are obtained using the fit 1/g^ Z -\-aq (left) and 
as an average over an interval in (p' (middle) as described in the text. For the final estimates, shown in the last 3 columns, we 
take the average of the two methods and half the difference (rounded up) for the errors. 



0.85^.... 

0 1 2 3 4 5 

q (GeV) 



q (GeV) 


FIG. 3. Data for Zt (upper) and ZtIZv (lower) after trans¬ 
lation to the MS scheme at 2 GeV as a function of the lattice 
momentum q . The lattice calculation was done on five ensem¬ 
bles in the RI-sMOM scheme. The a — 0.12 fm (a = 0.09 fm) 
fit are to the combination of al2m310 and al2m220 {al2m310 
and al2m220) data as there is no detectable dependence on 
the quark mass. The a — 0.06 fm fit is to the a06m310 en¬ 
semble data. The data were fit using the ansatz c/q^ + Z + aq 
and the dot-dashed, dotted and dashed lines show the ex¬ 
trapolation Z + aq for the a — 0.06, 0.09 and 0.12 fm data. 


offer two possible explanations. One, the large variation 
observed in the bin mean indicates that the ensembles 
of 0(1000) configurations (spanning a total of 5000-6000 
molecular dynamics evolution steps in the generation of 
thermalized HISQ lattices we have used) have not cov¬ 
ered enough phase space and bin errors are consequently 


underestimated. The other explanation is that, since we 
used the same four or eight source positions on all con¬ 
figurations in an ensemble, the data for fixed source po¬ 
sition is more correlated. Our ongoing tests confirm that 
using random but well-separated source positions on each 
configuration is a better strategy. Finally, based on the 
convergence of the bin distributions to a gaussian on in¬ 
creasing the bin size to the full sample and the lack of 
evidence of long tails, makes us confident that the final 
error estimates are reliable. 

Our overall conclusion about statistics is that while 
0 ( 10 , 000 ) measurements on these ensembles of 0 ( 1000 ) 
configurations are sufficient for extracting the tensor 
charge with few percent uncertainty, one will need a fac¬ 
tor of ten or more in statistics for obtaining the scalar 
charge with similar accuracy. This goal is currently being 
pursued using the AMA method discussed in Section [Tv| 

Lastly, we performed both correlated and uncorrelated 
fits to the nucleon two-point function data. In all cases in 
which the correlated fits were stable under changes in the 
fit ranges and had reasonable the two fits gave over¬ 
lapping estimates. Since correlated fits did not work in 
all cases, all statistical errors in the two- and three-point 
correlation functions were, thereafter, calculated using 
a single elimination jackknife method with uncorrelated 
fits performed on each jackknife sample. 


III. CONTRIBUTION OF THE CONNECTED 
DIAGRAM 


Estimates of the bare and renormalized charges on 
the nine ensembles at different lattice spacings, light 
quark masses and lattice volumes are given in Tables |VI| 
and IVIII 

To extrapolate these estimates to the physical point, 
he., the continuum limit (a —)■ 0 ), the physical pion mass 
{Mj^o = 135 MeV) and the infinite volume limit {L —>■ 
(xi), we explored the four parameter ansatz 


9t — 


1 + 


Mi 


(VF,)- 


ri^ 


H-C 2 G “h C 3 (/x) H- C 4 e 


-M^L 


( 12 ) 


where we have included the leading chiral logarithms m- 
The loop functions /'(/t/Mtt) for the two isospin channels 
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are: 






(2 + 4:g\) log + 2 + g\ 


2.72+ 6.38log ^ 


(2 + 8g\) log + 2 + 3g\ 


1.72+ 3.75 log 


(13) 


(14) 


where we use g, = Mp = 770 MeV for the renormaliza¬ 
tion scale and gA = 1.276. The extrapolation ansatz is 
taken to be linear in a because the discretization errors 
in the clover-on-HISQ formalism with unimproved opera¬ 
tors start at 0{a). Similarly, we have kept only the lead¬ 
ing finite volume correction term, . In Fig. we 


compare the fit obtained using Eq.(12| with that using 


the simpler isospin independent four parameter ansatz 
without the chiral logarithm: 


fits to the isovector — g^ and the connected part of 
the isoscalar g^ + g^ data using Eq. (15), are shown in 
Fig.|6l Our final estimates are 


gl^-'^{con) = 1.020(76), 

g^+'^icon) = 0.541(62). (17) 


with a x^/dof = 0.4 and 0.2 respectively. 


IV. CONTRIBUTION OF THE 
DISCONNECTED DIAGRAM 

In Section III Bl we showed that to estimate the discon¬ 
nected contribution, we need to calculate two quantities 
at zero-momentum—the nucleon two-point function and 
the contraction of the bilinear fermion operator into a 
quark loop—and measure their correlation. These two 
calculations are described below. 


gT{a,MT,,L) = Cl + C2a + csM^ + C4e . (15) 


A. Two-point Function 


Both fits have reasonable y^/dof and the estimates at 
the physical point are consistent. The fit including the 
chiral logarithm would naively indicate that gT should 
decrease in value with increasing for > 300 MeV. 
Such a behavior is not seen in the global data shown 
in Fig. We conclude that the large curvature due to 
the chiral logarithm seen in Fig. is most likely due to 
the number and accuracy of the data and of keeping just 
the leading chiral correction. Also, the error estimate 
from the fit using the simpler ansatz given in Eq.( |15[ ) is 
more conservative and covers the full range of both fits. 
We, therefore, use Eq. (15) for all further analysis in this 
paper. 

The results of the fits using Eq. © and the extrapo¬ 
lated value are shown in Fig. separately for operator 
insertion on the u and d quarks in the nucleon. We find 
that the g^ contribution is larger and essentially flat in all 
three variables (lattice spacing, pion mass and volume), 
while the gi^ connected contribution is much smaller and 
shows a slightly larger relative spread. The spread in 
g^ on the a = 0.06 fm lattices is an example of the un¬ 
expectedly large statistical fluctuations we mentioned in 
section El that will require higher statistics to resolve. 
The final renormalized extrapolated values for the proton 
charges are 


The high statistics calculation of the two-point func¬ 
tion with smeared sources was redone using the all-mode- 
averaging (AMA) technique |T7j because quark propa¬ 
gators from the earlier connected three-point function 
study were too expensive to store. To implement AMA, 
we again choose four different source time slices sepa¬ 
rated by Lt/ 4 on each configuration. On each of these 
time slices we calculate the two-point correlator by plac¬ 
ing Vlp = 15 low precision (LP) sources, for a total of 
4 X 15 = 60 sources per configuration. This estimate is a 
priori biased due to the LP calculation. In addition, on 
each of these four time slices we place one high precision 
(HP) source, be., IVhp = 4 such sources per configura¬ 
tion, from which we calculate a LP and a HP correlator. 
These four HP and LP correlators are used to correct the 
bias in the 60 LP estimates, be., on each configuration, 
the two-point function is given by 

-| 

-P(t,to) = ^ E 

i—1 

+ E [‘^HP ^0, xf P) - (t; to, xf P) , 

i—1 

(18) 


gi^{con) = 0.774(65), 

gi^icon) = -0.233(25). (16) 

The y^/dof is 0.1 and 1.6 for gij. and g^, respectively, 
with dof =5. In performing the fits, we assume that 
the error in each data point has a Gaussian distribution 
even though the quoted Ict error is a combination of the 
statistical error and the systematic error coming from 
the calculation of the renormalization factor Zt- The 


where and are the two-point correlation func¬ 
tion calculated in LP and HP, respectively, and and 
are the two kinds of source positions. 

The basic idea of AMA is that, in the low-precision 
evaluation, the LP average (first term in Eq. @) is bi¬ 
ased, and this bias depends predominately on low modes 
of the Dirac operator that are independent of the source 
position and can be corrected by the second term. Thus, 
we get an unbiased estimate from 60 LP source points for 
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FIG. 4. Comparision of the simultaneous fits versus a, and Mt^L to the isovector charge, data using Eq. ( |12[ | (top) 

with the simpler version without the chiral logarithms given in Eq. ( |15| l (bottom). The data symbols are defined in Table 
The fit is given by the red line and the physical value after extrapolation to the continuum limit (a —>■ 0), physical pion mass 
{Mtt —>■ M^o^°) and infinite volume (L —^ oo) is marked by a red star. The error band is shown as a function of each variable 
holding the other two at their physical value. The data are shown projected on to each of the three planes. 


ID 

con.ii 

con.d 

5t 

con,14 —d 

5t 

con,tt + d 

Ot 

disc,Z 

disc.s 

5t 

con,14 —d 

9v 

al2m310 

al2ni220S 

al2m220 

al2m220L 

0.875(18) 

0.873(26) 

0.888(22) 

0.859(18) 

-0.2208(93) 

-0.212(17) 

-0.222(12) 

-0.198(10) 

1.096(21) 

1.086(26) 

1.111(24) 

1.058(19) 

0.655(20) 

0.661(36) 

0.665(26) 

0.662(21) 

-0.0124(23) 

-0.0038(41) 

-0.0041(20) 

-0.0010(28) 

1.069(9) 

1.059(12) 

1.074(11) 

1.065(7) 

a09m310 

a09m220 

a09ml30 

0.829(16) 

0.820(16) 

0.779(33) 

-0.2025(80) 

-0.2120(79) 

-0.214(18) 

1.031(19) 

1.033(17) 

0.993(33) 

0.626(18) 

0.608(18) 

0.565(42) 

-0.0050(21) 

-0.0005(21) 

-0.0021(53) 

1.056(8) 

1.050(9) 

1.029(16) 

a06m310 

a06m220 

0.778(18) 

0.759(43) 

-0.1898(86) 

-0.241(19) 

0.969(20) 

1.002(46) 

0.588(21) 

0.519(48) 

-.0035(62) 

-0.0005(52) 

1.040(8) 

0.993(18) 


TABLE VI. The bare connected (<?“") and disconnected contributions to the tensor charges of the proton on the nine 

ensembles. Dashes indicate that those ensembles have not been simulated. The isovector vector charge gY~‘^ is used to construct 
ratios for noise reduction as described in the text. 


the computational cost of (60 + 4) LP and 4 HP calcu¬ 
lations. In our current implementation, 15 LP measure¬ 
ments cost the same as one HP when using the multigrid 
algorithm for inverting the Dirac matrix |22j . (On the 
aOSmS 10 ensemble we used (92-1-4) LP and 4 HP sources 
and the errors decreased by a factor of ^ 1.2 compared 
to (60-1-4) LP sources.) Comparing the errors in the es¬ 
timates for masses given in Table |IV[ we find that the 
AMA errors are a factor of 2 — 4 times smaller than those 
from the connected study (all HP measurements). Since 
this improvement is based on comparing 120 LP (we ef¬ 
fectively doubled the LP statistics by analyzing both the 
forward and backward propagation of the nucleon) versus 
8 HP measurements, we conclude that the correlations 
between the LP measurements are small. 

If we assume that the variance of both the LP and HP 


measurements is the same and given by a and the cor¬ 
relation between the HP and LP measurements from the 
Ahp points, C = crJfpLp/^'^^ i® small, then the statistical 
error in Eq. (181 is given by m 


The second term in the square root becomes smaller as 
the LP estimate approaches the HP estimate and the cor¬ 
relation factor C —1. By controlling A^lp, A^hp and C, 
we can minimize the total error for a fixed computational 
cost. 

To speed up the AMA method we exploit the fact 
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ID 

_con,ii 

9t 

coxi.,d 

9t 

_con,u — d 

9t 

con,ii-|-d 

9t 

_discJ 

9t 

_disc,s 

9t 

al2m310 

0.852(37) 

-0.215(12) 

1.066(46) 

0.637(31) 

-0.0121(23) 

-0.0040(19) 

al2m220S 

0.857(43) 

-0.209(19) 

1.066(50) 

0.649(44) 

— 

— 

al2m220 

0.860(40) 

-0.215(15) 

1.075(48) 

0.644(36) 

-0.0037(40) 

-0.0010(27) 

al2m220L 

0.840(37) 

-0.194(12) 

1.033(45) 

0.647(33) 

— 

— 

a09m310 

0.840(28) 

-0.2051(98) 

1.045(34) 

0.634(25) 

-0.0050(22) 

-0.0005(21) 

a09m220 

0.836(28) 

-0.216(10) 

1.053(34) 

0.619(25) 

— 

-0.0021(54) 

a09ml30 

0.809(40) 

-0.222(20) 

1.032(44) 

0.587(45) 

— 

— 

a06m310 

0.815(29) 

-0.199(10) 

1.015(34) 

0.617(27) 

-0.0037(65) 

-0.0005(55) 

a06m220 

0.833(52) 

-0.264(22) 

1.099(59) 

0.569(55) 

— 

— 


TABLE VII. The renormalized connected (gr) and disconnected contributions to the tensor charges of the proton on the 

nine ensembles. The errors are obtained by adding in quadratures the statistical errors given in Table [Vl| in the bare charges 
to the errors in the renormalization constants given in Table [y] 
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FIG. 5. Simultaneous extrapolation to the physical point (a —>■ 0, Mtt —>■ and L —>■ oo) using Eq. | |15[ ), of the connected 

contributions to the flavor diagonal nucleon (proton) tensor charges, (upper) and g^ (lower), renormalized in the MS scheme 
at 2 GeV. The physical values given by the fit are marked by a red star. The rest is the same as in Fig. 


that the same Dirac matrix is inverted multiple time^ 
on each configuration. It is, therefore, efficient to pre¬ 
condition the matrix by deflating the low-eigenmodes. 
We implement such improvement using the multigrid 
solver [22l [23] which has deflation built in. To obtain 
the LP estimate of the two-point function, we truncate 
the multigrid solver using a low-accuracy stopping cri¬ 
terion: the ratio (rup = | residue |lp/| source |) is chosen 
to be 10“^ for all the ensembles. Our final analysis of 
the masses, amplitudes and matrix elements, however, 
shows that this stopping criteria was overly conservative 
as the bias correction term is negligible compared to the 
statistical errors. 


^ The number of inversions of the Dirac matrix per configuration 
are 12 X (A^lp -f Vhp) 


B. Disconnected Quark Loop 


For the evaluation of the quark loop term 
Tr[M“^(T, x; r, x)r], we adopt the stochastic 
method accelerated with a combination of the truncated 
solver method (TSM) [TS] [TB] , the hopping parameter 
expansion (HPE) [24l [23 and the dilution technique 
[26ll28j . To obtain a stochastic estimate of the quark 
loops, consider a set of random complex noise vectors \r]i) 
for f = 1, 2, 3, • • • , N, having color, spin and spacetime 
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FIG. 6. Simultaneous extrapolation to the physical point (a —>■ 0, —>■ and L —>■ oo) using Eq. (15 I of the connected 

contributions to the isovector (upper) and isoscalar (lower) nucleon (proton) tensor charges renormalized in the MS 
scheme at 2 GeV. The overlay in the middle upper figure, with the dashed line within the thin grey band, is the fit to the data 
versus assuming no dependence on the other two variables. Rest is the same as in Fig. 


components with the following properties: 


by 


N 


Vn 


i—l ^ ^ 


( 20 ) 

( 21 ) 


We choose complex Gaussian noise vectors, i.e., we fill all 
the spin, color and spacetime components of the vector 
with (vr + iri)/\/2, where and are Gaussian random 
numbers, because they give marginally smaller statistical 
error than Zjy random noise when combined with the 
HPE. 

These random vectors are used as sources for the in¬ 
version of the Dirac matrix. Then, from the solutions 
I Si) of the Dirac equation, 


M\si) = \r]i) , 

the inverse of the Dirac matrix is given by 


( 22 ) 


M 


-1 


1 ^ / 1 ^ \ 
= ]v ^ ~ iv ^ ) 

2^1 \ i^l ) 




(23) 

(24) 


The stochastic estimate of the zero-momentum insertion 
of the operator contracted into a quark loop is then given 


N 


^Tr[M i(T,x;T,x)r] « — ^(7?,|^r|si)^, (25) 

X i—l 

where \rj)^ is a vector whose t = t components are filled 
with random numbers, and the entries on other time 
slices are set to zero. 

In the estimation of the inverse of the Dirac matrix by 
using random sources, Eq. ( [24| , one can use the mixed- 
precision technique called the truncated solver method 
(TSM) [13 [T2]. The idea of the TSM is the same as the 
AMA used in the evaluation of the two-point function. 
Consider two kinds of solution vectors of Eq. (22) for a 
given random source \rji) with different precision: |si)LP 
and |si)HP, where |si)LP is the low precision computa¬ 
tionally cheap estimate of the solution, while |si)HP is 
the high precision solution. The low precision estimate, 
|si)LP, was obtained by truncating the multigrid inverter 
at rpp = 5 X 10“^. The bias observed with this choice of 
rpp will be discussed later in this section. 

By using the LP and HP solutions, the unbiased esti¬ 
mator of M~^ is again given by 

A^lp+AThp 

+ ^ - X (ki)HP - |Si)Lp) (??i| • (26) 

i=ArLP-|-l 

The first term in the right-hand-side (r.h.s) is the LP es¬ 
timate of M~^ while the second term in the r.h.s corrects 


= 

-/Vlp 
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the bias. As described in the case of the two-point func¬ 
tion estimation, the total statistical error of scales 
as Eq. (19). In other words, there are again two sources 
of statistical error in one is the LP estimate that 


scales as 1 /A^lp ; and the other is the correction term 
that scales as -^/I/A^hp- The size of the statistical error 
in the correction term is determined by the correlation 
between |si)HP and |si)Lp. 

In the TSM, there are three parameters we can tune 
to minimize the statistical error for a given computation 
cost: Alp, Alp/Ahp and the LP stopping criteria tlp. 
Note that once Alp/Ahp and clp are determined, the 
total error scales as -y/l/ALp. Hence Alp determines the 
size of the error, and Alp/Ahp and rpp determine the 
efficiency of the estimator in terms of the computational 
cost. To maximize the efficiency, we tune the App/Aup 
and rpp so that the size of the error from the correction 
term is much smaller as it minimizes the computation 
time. In th is study, we use App/Anp = 30 or 50 (See 
Table. |VIII ) and rpp ^ 5 x 10“^. With this accuracy, 
we find that the bias correction term is ~ 10 % of the 
final estimate of gip{disc) and about half of the statistical 
error. 

We improve the TSM by using the hopping parameter 
expansion (HPE) [24l |25] as a preconditioner to reduce 
the statistical noise. In the HPE one writes the clover 
Dirac matrix as 


M = — (1 - kD) , 


(27) 


where k is the hopping parameter. The inverse can then 
be written as 


- n— 1 


2=1 




(28) 


By taking n = 2, the disconnected quark loop is given by 

Tr [M-^T] = Tr [(2 k 1 -h 2k^D + P] . (29) 

Here, the first two terms of the r.h.s do not contribute to 
the nucleon tensor charge because TrP = Tr(r 7 ^) = 0. 
As a result, the only nontrivial term that we need to 
calculate is Tt Because the two leading 

terms, which would otherwise contribute only to the 
noise, are removed from the stochastic estimation, HPE 
works as an error reduction technique. Tests using the 
al2m310 ensemble show that the statistical error of the 
disconnected contribution to the tensor charge is reduced 
by a factor of about 2.5 with HPE. 


As shown in Eq. (23), the noise in the stochastic es¬ 
timation for M~^ is proportional to whose mag¬ 

nitude decreases exponentially as the spacetime distance 
between source and sink increases. Hence it is possible 
to reduce the statistical noise by placing noise sources 
only on part of the whole time slice, choosing maximally 
separated points, and fill the other points on the time 
slice with zero. This procedure divides the time slice 


into m subspaces, and the answer for the full time slice is 
obtained by combining results of the m subspaces. The 
computational cost increases by a factor of m because 
Dirac inversions are needed for each noise source vector 
defining a subspace. Hence this technique is useful when 
the reduction in noise wins over the increase in computa¬ 
tional cost. This is called the time dilution method ng- 
[25j. Unfortunately, we find that the increase in com¬ 
putational cost is equal to or larger than the gain from 
the reduction of statistical noise for the nucleon charges. 
Hence we place random sources on all points of a time 
slice and for each time slice that we want to evaluate the 
operator on. 

There is one more symmetry that can be used for noise 
reduction: 75 -hermiticity of clover Dirac operator, Afl = 
75 M 75 . Because of this symmetry, the quark loop for 
tensor channel should be pure imaginary, and the nucleon 
two-point function is real. Hence we set the real part of 
the quark loop to zero when constructing the correlation 
function and averaging over the configurations in Eq. ([^ . 

To increase the statistics, we average over the three 
possible combinations of 7 ^ 7 ^ and forward/backward 
propagators. The final values of tsep investigated, the dis¬ 
placement r with respect to the source time slice of the 
two-point correlator on which the operator was inserted, 
the statistics and the number of random noise sources 
used on each configuration are given in Table |VHI| 


C. Results for the disconnected contributions 

The calculation of the disconnected diagram is com¬ 
putationally expensive so it has been done on only four 
ensembles: al2m310, al2m220, a09m310 and a06m310. 
These four ensembles provide an understanding of the 
discretization errors and of the behavior as a function 
of the quark mass. To get the full contribution of the 
quark EDM to the nucleon EDM, we also need to evalu¬ 
ate the disconnected diagram with a strange quark loop. 
Since the calculation with the strange quark are compu¬ 
tationally cheaper, we have also analyzed a fifth ensem¬ 
ble, a09m220, for that estimate. 

We use the fit ansatz given in Eq. ( [II| ), i.e. the same 
two-state-fit method used to fit the data for the connected 
three-point diagrams, to extract the ground state results 
for the disconnected contribution. The data and the re¬ 
sults of the fit for the light and strange quark loop on 
the al2m310 ensemble are shown in Figs. and re¬ 
spectively. We find significant contribution from excited 
states only on the al2m310 ensemble for light quark dis¬ 
connected diagram — it is large for tgep = 8 , but by 
tsep = 12 the data agree with the final extrapolated value. 
The peculiar pattern seen in the a06m310 ensemble is 
most likely due to the small number (100 [ 200 ]) of con¬ 
figurations analyzed as given in Table |VHl| 

Including the disconnected diagrams also requires cal¬ 
culating their contribution to the renormalization con¬ 
stants in the RI-sMOM scheme. We have not done this 
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FIG. 7. Fits, using Eq. (Ill, to isolate the excited state contribution in the light quark disconnected diagram, are shown 

for the four ensembles analyzed. The solid black line and the grey band are the ground state estimate and error. The data and 


results of the fit for different taep are also shown. 
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FIG. 8. Fits, using Eq. (Ill, to isolate the excited state contribution in the strange quark disconnected diagram, 
shown for the five ensemh 


rles analyzed. The solid black line and the grey band are the ground state estimate and error. The 


due to the poor signal in disconnected diagrams and use 
the same renormalization factor as calculated for the con¬ 
nected diagrams. In perturbation theory, the discon¬ 
nected diagrams come in at higher order, so their contri¬ 
butions are expected to be small. Furthermore, the dis¬ 
connected quark loop contributions themselves are very 
small for the nucleon tensor charges, so we expect the 
impact of the small difference in the renormalization fac¬ 
tor due to neglecting the disconnected piece in Zt will 
change the final estimate by much less than the statistical 


error quoted in Table VI The hnal renormalized results, 
with this caveat, are given in Table |VlT] 


There are two ways in which we can study the quark 
mass dependence of the disconnected contribution. First, 
by comparing the strange with light quark loop contri¬ 
butions we note that the estimates on all four ensembles 
increase as the quark mass is decreased. The second is to 
compare the estimates on the al 2m310 a.nd al2m220 en¬ 
sembles. Unfortunately, the statistical errors in the latter 
are too large to draw a conclusion, even though we used 
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the largest number of random sources per configuration 
for this study. Our conclusion is that a higher statistics 
study is needed to quantify the quark mass dependence 
and reduce the overall error in the disconnected contri¬ 
bution so that a reliable continuum extrapolation can be 
made. 

The authors in Ref. [H] found that the disconnected 
contribution to the nucleon tensor charge is consistent 
with zero on a Nf = 2 -|- 1 -|- 1 twisted mass fermion en¬ 
semble at a = 0.082(4) fm and Mjr = 370 MeV. While a 
direct comparison with our results would be meaningful 
only after both results have been extrapolated to the con¬ 
tinuum and physical pion mass limit, we note that our 
estimates are also consistent with zero for all ensembles 
with the strange quark loop, and in two of the four cases 
of light quark loops. 

Given that the estimates of the disconnected contribu¬ 
tion with light quark loops are small compared to con¬ 
nected part, have large errors, and have been obtained 
on only four ensembles, we do not include them in esti¬ 
mates of the isoscalar charges Instead, we take the 

largest value 0.0121 on the alEmSlO ensemble and use it 
as an estimate of the systematic error associated with 
neglecting the disconnected piece. This error is added 
in quadrature to the overall error in the connected esti¬ 
mate. The disconnected contribution with strange quark 
loops is even smaller but we keep it since it does not have 
a connected piece and we can perform a reasonable ex¬ 
trapolation in the lattice spacing and the quark mass as 
shown in Fig. and get 

^Misc ^ 0.008(9), (30) 

with a x^/dof = 0.29 for dof = 2. Bounding is impor¬ 
tant for the analysis of the neutron EDM, especially if the 
chirality flip is controlled by the Higgs Yukawa coupling. 
In those BSM scenarios, the contribution of would 
be enhanced by the ratio of quark masses rns/rnu,d (be., 
proportional to the coupling of a “Higgs” to quarks), rel¬ 
ative to g^ and g^. Using these estimates, the analysis 
of the contribution of the quark EDMs to the neutron 
EDM in presented in Section |VB| 


V. NUCLEON TENSOR CHARGES AND 
QUARK ELECTRIC DIPOLE MOMENT 


In the previous Sections [Ill| and [TV C[ we discussed the 
calculation of the connected and disconnected diagrams 
to the nucleon tensor charges. In this section we present 
our final results for the nucleon tensor charges and the 
constraints they put on the quark EDM couplings using 
the current bound on the neutron EDM. 


A. Nucleon Tensor Charge 

The isovector tensor charge needed to probe 

novel tensor interactions at the TeV scale in the helicity- 


flip part of neutron decays, does not get any contributions 
from the disconnected diagram in the iso-spin symmetric 
limit that we are working under. We consider the ex¬ 
traction of g^~‘^ reliable because all systematics are un¬ 
der control. In particular, we And (i) that the fit ansatz 
in Eq. (II) converges, indicating that the excited state 
contamination has been isolated, (ii) the data for the 
renormalization constant in the RI-sMOM scheme shows 
a window in for which the final estimates in the MS 
scheme at 2 GeV are constant within errors as discussed 
in Section HD Finally (iii), the estimates on the nine 
ensembles show little dependence on the lattice spacing, 
pion mass and lattice volume as shown in Fig. and dis¬ 
cussed in Section iHll 

Our final estimate given in Eq. ( |l7| ), g^~‘^ = 1.020(76), 
is in good agreement with other lattice calculations 
by the LHPC {Nf =2-1-1 HEX smeared clover ac¬ 
tion, domain wall action, and domain wall-on-asqtad ac¬ 
tions) [30] j RBG/UKQCD {Nf = 2-1-1 domain wall 
fermions [3T], ETMC {Nf = 2 -|- 1 -|- 1 twisted mass 
fermions) and the RQGD {Nf = 2 0{a) improved 

clover fermions) [33] as shown in Fig. A more de¬ 
tailed discussion of the systematics in these calculations 
and control over them using the FLAG quality criteria [6] 
is given in the Appendix [A| 

An analysis of the extrapolation to the physical quark 
mass has also been carried out by the LHPC [30] and 
RQGD [33] collaborations. They did not find signif¬ 
icant dependence on the lattice spacing and volume, 
so they extrapolate only in the quark mass using lin¬ 
ear/quadratic (LHPC) and linear (RQGD) fits in M"^. 
Their final estimates, g^'^ = 1.038(11)(12) (LHPC) and 
= 1.005(17) (29) (RQGD) are consistent with ours, 
but the size of our error is much larger. This is due to 
a combination of three factors in our calculation: (i) our 
determination of renormalization constants have larger 
uncertainty; (ii) errors in individual points are larger be¬ 
cause they are the estimates in the t, 


sep 


oo limit ob¬ 
tained by extrapolating the data with multiple Gep using 
a two state ansatz; and (iii) we extrapolate in all three 
variables using Eq. (|T^ , whereas LHPC and RQGD ex¬ 
trapolate only in M^~A fit to our data versus only , 
also shown in Fig gives a similarly accurate estimate 
g'^~'^ = 1.059(29) with a x^/dof = 0.3. 

A comparison between recent lattice QCD results for 
g^~^ and estimates derived from model ca lculat ions and 
experimental data are summarized in Fig. Ipp The lat¬ 
tice estimates show consistency and little sensitivity to 
the number of flavors, be., Nf = 2 or 2 -|- 1 or 2 -|- 1 -I- 1, 
included in the generation of gauge configurations. The 
errors in model and phenomenological estimates (integral 
over the longitudinal momentum fraction of the experi¬ 
mentally measured quark transversity distributions) are 
large. Only the Dyson-Schwinger estimate (DSET4) has 


^ A similar comparison is presented in Ref. m- 
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a (fm) (GeV^) 

FIG. 9. Fits using Eq. ( |15| l to obtain the result in the continuum limit (a —> 0) and at the physical pion mass (M.^- —>■ 
of the strange quark disconnected contribution. A finite volume study was not carried out for the disconnected contribution. 
Rest is the same as in Fig. 
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^sep 

/a 
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■p.Tdisc,l 
^ coni 

'^conf 

-^^LP 

A rdisCjS 

-'^LP 

Npp/Npip 

al2m310 

{8~ 

14} 

3,4,--- ,11 

1013 

1013 

5000 

1500 

30 

al2m220 

{8~ 

14} 

3,4,--- ,11 

958 

958 

11000 

4000 

30 

a09m310 

{10^ 

M6} 

6, 7, 8, 9 

1081 

1081 

4000 

2000 

30 

a09m220 

{10^ 

M6} 

5, 6, 7, 8,9 

— 

200 

10000 

8000 

50 

a06m310 

{16^ 

24} 

6,8,10,-,18 

100 

200 

10000 

5000 

50 


TABLE VIII. The parameters used in the study of the disconnected diagrams. The source-sink time separations analyzed (tsep), 
the time slices (r) on which the operator is inserted as explained in the text, the number of configurations analyzed (Ncont) and 
the number of random noise sources (Nlp) used on each configuration. Here {A ~ B} denotes the set of consecutive integers 
from A to B. The number Npp/Nnp gives the ratio of the number of low to high precision calculations done. The LP criteria 
for stopping the Dirac matric inversion was set to rpp = 0.005. For the associated two-point function calculation, we used 
AMA with 64 LP and 4 HP measurements on each conhguration and the results for the masses and amplitudes are given in 
Table HVl 


comparable errors and lies about da below the lattice 
QCD estimates. 

To summarize, even with a very conservative error es¬ 
timate, our result = 1.020(76), meets the target 

uncertainty of ~ 10% required to bound novel tensor in¬ 
teractions using measurements of the helicity flip part of 
the neutron decay distribution in experiments planning 
to reach 10“^ accuracy. Our goal for the future is to re¬ 
duce the error in gs, which currently is ~ 30% for the 
data sets presented in this work, to the same level. 


B. Quark Electric Dipole Moment 


The calculation of the connected contr ibution to the 
g!^ has been discussed in Section V A Estimates of 
the discon nected contribution were discussed in Sec¬ 
tion IV C Including the largest value (0.0121 obtained 
on the al2m310 ensemble) as a systematic error, our fi¬ 
nal results in the MS scheme at 2 GeV for the nucleon 
charges that get contributions from the disconnected di¬ 
agrams are: 


= 0.774(66), 
g^ = -0.233(28), 

g^+<^ = 0.541(67). (32) 


The quark EDM contributions to the neutron EDM, 
dm are given by 

dn = d-it T d-d gx + g^ (31) 

where the low-energy effective couplings dm da and dg 
encapsulate the new CP violating interactions at the TeV 
scale. The goal of the analysis, knowing the charges g^. 
and a bound on is to constrain the couplings dq and, 
in turn, BSM theories. 


Note that incorporating the disconnected contribution as 
a systematic error increases the errors marginally as can 
be seen by comparing estimates in Eq. (32) with those 
in Eq. (16). Results for the neutron tensor charges are 
obtained by using the iso-spin symmetry, i.e., by inter¬ 
changing the labels u ^ d. 

These final estimates are significantly smaller in mag¬ 
nitude than the quark model values, g^ = 4/3 and 
g^ = —1/3, but consistent with estimates derived from 
model calculations and experimental data summarized in 
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FIG. 10. A comparison between recent lattice QCD results 
for and estimates derived from model calculations and 
experimental data. The published lattice QCD results are 
from LHPCT2 |30], RBC/UKQCDTO [3l], RQCDT4 [35] 
and RBC’08 [^. Lattice estimates with reasonable con¬ 
trol over excited state contamination and extrapolation to 
the physical pion mass and the continuum limit are shown 
in green. Estimates from models and phenomenology are 
from Bacchetta’13 [38], Anselmino’13 [39], Kang’15 [40], Sum 
Rules’OO [H], DSE’14 [H]. 


FIG. 11. Bounds on the couplings du,d for the case = 0. 
Estimates used for and are given in Eq. (32|. 


develop the lattice methodology to include the contribu¬ 
tions of the quark chromo electric dipole moment opera¬ 
tor. 


Table jlXPj T he three phenomenological estimates Bac- 
chetta 13 |38j. Anselmino’13 [33], and Kang’15 [3D] give 
consistent but lower estimates for and g^ with g^~'^ ~ 
0.65. Similarly, taking the errors at face value, the 
Schwinger-Dyson estimate is ~ 4cr below the lattice QCD 
results. A recent reevaluation of the calculation of ten¬ 
sor charges using QCD sum rules with input from lattice 
QCD has been reported in [44l|45]. Their estimates in 
the MS scheme at 1 GeV are g^ = 0.79 and g^ = —0.20, 
each with « 50% uncertainty. Run to 2 GeV, these es¬ 
timates would decrease by « 10% in magnitude. These 
results are consistent with ours given in ( j^ but place 
less stringent constraints on the neutron EDM and BSM 
theories due to the larger uncertainty. 

Assuming that only the EDMs of the u, d, and s quarks 
contribute to the neutron EDM via Eq. (31) and the val¬ 
ues of g^’ are given by Eqs. (321 and (301, one can 
put bounds on the du,d,s- Using the current estimate 
Idffl < 2.9 X 10“^® e cm (90% CL) [33], 1-sigma slab pri¬ 
ors for g^ and g^ given in Eq. ( |^ , and assuming g^ — 0^ 
we obtain the 90% confidence interval bounds for and 
dd shown in Fig. El Note that dg is not constrained since 
g^ is consistent with zero. 


Using these estimates of we have analyzed the 

consequences on split SUSY models, in which the quark 
EDM is the leading contribution in [5]. Our goal for the 
future is to improve the estimates presented here and 


® The effect of the choice of scale fi is illustrated by converting the 
results in Ref. m from = 1 GeV to /i = 2 GeV using two-loop 
running with the anomalous dimensions taken from Refs. |20II43| . 


VI. CONCLUSIONS 


We have presented a high statistics study of the isovec¬ 
tor and isoscalar tensor charges of the nucleon using 
clover-on-HISQ lattice QCD. We calculate both the con¬ 
nected and disconnected diagrams contributing to these 
charges. The analysis of nine ensembles covering the 
range 0.12 — 0.06 fm in lattice spacing, M,r = 130 — 
320 MeV in pion mass, and M-j^L = 3.2 — 5.4 in lattice 
volume allowed us to control the various sources of sys¬ 
tematic errors. We show that keeping one excited state 
in the analysis allows us to isolate and mitigate excited 
state contamination. The renormalized estimates of the 
various tensor charges show small dependence on the lat¬ 
tice volume, lattice spacing and the light quark mass. 
These results can, therefore, be extrapolated reliably to 
the physical point. 

Our final estimate for the tensor charge g^~‘^ = 
1.020(76) is in good agreement with previously reported 
estimates. The signal in the calculation of the discon¬ 
nected diagrams is weak in spite of using state-of-the-art 
error reduction techniques. The value is small and we 
bound its contribution to light quark charges g^ and g^. 
The signal for strange disconnected loop is even smaller, 
however in this case we are able to extrapolate the results 
to the continuum limit and find g^ = 0.008(9). Using 
these estimates and the current bound on the neutron 
electric dipole moment, we carry out a first lattice QCD 
analysis of the constraints on the strengths of the up, 
down and strange quark electric dipole moments. The 
impact of these constraints on the viability of split SUSY 
models, in which the quark EDM is the leading contri¬ 
bution to the neutron EDM, is carried out in [5]. 
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9t 

Qt 

5t 

g 

This study 

-0.23(3) 

0.77(7) 

0.008(9) 

2 GeV 

Quark model 
QCD Sum Rules fiT] 
Dyson-Schwinger [32] 
Transversity 1 [3S] 
Transversity 1 [3H] 
Transversity 2 [3D 
Transversity 3 [30] 

-1/3 
-0.35(17) 
-0.11(2) 
-0.18(33) 
-0.16(30) 
-0.25(20) 
—0 22+“'^"^ 

4/3 

1.4(7) 

0.55(8) 

0.57(21) 

0.51(19) 

0.39(15) 

0.39t°:?I 

- 

7 

2 GeV 

~ 1 GeV 

2 GeV 

~ 1 GeV 

3.2 GeV 


TABLE IX. A comparison of our lattice estimates of and of the proton with those from different models and phenomenol¬ 
ogy. The “Transversity 1” estimate is given both at the original scale at which it was evaluated (~ 1 GeV) and after running 
to 2 GeV to show the magnitude of the scaling effect. The symbol “?” in the last column indicates that the scale at which the 
calculation is done is undetermined. 
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Appendix A: Systematics in the calculation of the 
isovector nucleon tensor charge 

In Table we give a summary, in the FLAG for¬ 
mat m, of the level of control over various systematics 
in the calculation of the isovector tensor charge of the nu¬ 
cleon using simulations of lattice QCD with Nf = 2, 2-1-1 
and 2 -I- 1 -I- 1 flavors. Note that a community wide con¬ 
sensus on applying the FLAG criteria to matrix elements 
within nucleon states does not yet exist. By performing 
this analysis, we wish to emphasize that the agreement 
between various calculations of g^~‘^ has reached a level 
of precision that calls for a FLAG like analysis. 

The systematics covered by the FLAG criteria are also 
encountered in the calculation of matrix elements within 
baryon states. We, therefore, follow the same quality 
criteria for the publication status, chiral extrapolation, 
finite volume effects, and renormalization as defined by 
FLAG [6] and define an additional criterion, excited state 
contamination, that is relevant to the calculations of ma¬ 
trix elements within nucleon states. For the criterion 
“continuum extrapolation” we relax the requirement of 
an extrapolation provided the data meet the rest of the 
requirements: do not warrant an extrapolation, and a 
reasonable estimate of the uncertainty is provided. We 
also do not require that the action and the operators are 
0{a) improved. 

• Publication status: 

A published or plain update of published results 

P preprint 

C conference contribution 

• Chiral extrapolation: 

★ A7,r,min < 200 MeV 

O 200 ’MeV < < 400 MeV 

■ 400 MeV < 
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Continuum extrapolation: 

★ 3 or more lattice spacings, at least 2 points 
below 0.1 fm 

1.2 

O 2 or more lattice spacings, at least 1 point 
below 0.1 fm 

1.1 

■ otherwise 

■a 


it^l.O 


bjo 

Finite-volume effects: 

A^ 7 i-,min.b > 4 or at least 3 volumes 

0.9 

O > 3 and at least 2 volumes 

■ otherwise 

0.8 


1—'—'—'—'—I—'—'—'—'—r 



_ i . PNDME 4f clover/HISQ T LHPC 3f clover ♦ETMC2fTMF J 

♦ ETMC 4f IMF T LHPC 3f DWF ► RQCD 2f clover . 

T LHPC 3f DWE/asqtad 
■ RBC/UKQCD 3f DWF 

j_^_,_,_,_I_,_,_,_,_I_,_,_,_^_: 

0 0.05 0.10 0.15 


• Renormalization: 

★ nonperturbative 

O 1-loop perturbation theory or higher with a 
reasonable estimate of truncation errors 

■ otherwise 

• Excited State: 

tsep, max > 1.5 fm Or at least 3 source-sink 
separations, tgepj investigated at each lattice 
spacing and at each 

O At least 2 source-sink separations with 
1.2 fm < tsep, max < 1.5 fm at at least one M^r at 
each lattice spacing. 

■ otherwise 


Plots of the data, summarized in Tabl e pQ as a function 
of a, and are shown in Fig. |12| One observes 
very little sensitivity to these three variables and on the 
number of fermion flavors or the lattice action used. 


a(fm) 


1.2 



Q _ • PNDME 4f clover/HISQ t LHPC 3f clover ♦ ETMC 2f TMF J 

^ ^ . ♦ ETMC 4f TMF t LHPC 3f DWF ► RQCD 2f clover . 

! T LHPC 3f DWF/asqtad 

I ■ RBC/UKQCD 3f DWF 

0.8 ^ . i .. I .... I .... I .... 1 

0 0.05 0.10 0.15 0.20 


Ml (GeV^) 


1.2 

1.1 


13 

it^l.O 

bjo 


0.9 


0.8 


-i—i—I—^—r 



hJh ^ 1 




• PNDME 4f clover/HISQ 
► ETMC 4f TMF 


LHPC 3f clover ♦ ETMC 2f TMF 

LHPC 3f DWF ► RQCD 2f clover , 

LHPC 3f DWF/asqtad 
RBC/UKQCD 3f DWF 

I ■ ■ ■ I ■ ■ ■ 


10 


M„L 


FIG. 12. Estimates of from lattice QCD for Nf = 2, 
2-1-1 and 2 -|- 1 - 1 - 1 flavors from the PNDMET5 (this work), 
ETMCT5 [32H34] LHPCT2 [30l, RBC/UKQCDTO [H], and 
RQCDT4 collaborations. These data show little sensitiv¬ 
ity to a (top), (middle), (bottom) and on whether 

the strange and charm quarks are included in the generation 
of the lattice ensembles or on the lattice action used. The 
vertical dashed line in the middle panel marks the physical 
pion mass M,r = 135 MeV. 
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Collaboration 


Ref. 














»9 








®S 






Nf 


















<eN>- 


9T 


PNDME’15 

This work 

p 

2+1+1 

★ 

★ 

★ 

★ 

★ 

1.020(76; 

a 

ETMC’15 

m 

c 

2+1+1 

■ 

■ 

★ 

★ 

★ 

1.053(2i; 

n 

LHPC’12 

m 

A 

2+1 

★ 

o 

★ 

o 

★ 

1.038(11' 

uTI 

RBC/UKQCD’IO 

m 

A 

2+1 

O 

■ 

★ 

■ 

★ 

0.9(2fl 


RQCD’14 

ESI 

A 

2 

★ 

★ 

★ 

o 

★ 

1.005(17) (29f 

ETMC’15 

m 

P 

2 

★ 

■ 

■ 

★ 

★ 

1.027(^ 

!l 

RBC’08 


A 

2 

■ 

■ 

★ 

■ 

★ 

0.93(6fl 



^ This estimate is obtained from a simultaneous fit versus a, M^, and defined in Eq. l |15| l using data on nine clover-on-HISQ 

ensembles given in Table IVIII 

^ The quoted estimate m is from a single Mtt = 373 MeV, a = 0.082 fm and Nf = 2 + 1 + 1 maximally twisted mass ensemble. Three 
values of tgep ~ 1? 1-15, and 1.3 fm are analyzed for handling excited state contamination. We quote their result from the two-state fit. 
A second low statistics study on an ensemble with Mtt = 213 MeV and a = 0.064 fm gave a consistent estimate. 

^ The central value is from a two parameter chiral fit to just the coarse Wilson ensembles data. This agrees with a three parameter 
chiral fit to data from three different lattice actions simulated at different lattice spacings and with different volumes. Uncertainty due 
to extrapolation in the lattice spacing a and the finite volume controlled by Mt^L is expected to be small. 

^ Result is based on simulations at one lattice spacing 1/a = 1.73 GeV using domain wall fermions. The statistics for the ensembles 
corresponding to the four pion masses simulated, Mtt = 329, 416, 550, 668 MeV, were 3728, 1424, 392, 424 measurements, 
respectively. A single tsep = 1-39 fm was used. 

® The result of this clover-on-clover study is obtained using a fit linear in keeping data with <0.1 GeV^ only. Data do not show 
significant dependence on lattice spacing or lattice volume. Excited state study is done on three of the eleven ensembles. Most of the 
data are with £sep ~ 1 fm. The second error is an estimate of the discretization errors assuming they are O(a^) since 0(a) improved 
operators with 1-loop estimates for the improvement coefficients are used in calculations done on a = 0.081, 0.071 and 0.06 fm lattices. 
Preliminary estimates presented by the QCDSF collaboration |48| are superseded by this publication m- 

^ Result from a single ensemble of maximally twisted mass fermions with a clover term at a = 0.093(1) fm, Mtt = 131 MeV and 
M-rrL ^ 3. To control excited state contamination, three values of tsep ~ 0.94, 1.1 and 1.3 fm are analyzed. We quote their value from 
the £sep ~ 1.3 fm analysis. 

s Results based on one lattice spacing 1/a = 1.7 GeV with the DBW2 domain wall action, three values of quark masses with 
Mtt = 493, 607, 695 MeV, and 0(500) measurements. Only one tsep = 10 (1.14 fm) was simulated except at the lightest mass where 
tsep = 12 data was generated but used only as a consistency check as it has large errors. 


TABLE X. A summary of the control over various sources of systematic errors in lattice QCD calculations of the isovector 
tensor charge using the FLAG quality criteria [6] reproduced in this Appendix. Results from all collaborations quoted in 
this table have used nonperturbative methods for calculating the renormalization constants. 
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